Acoustics in one dimension¶
This notebook demonstrates running the Clawpack Fortran code and plotting results from a Jupyter notebook, and illustration of the behavior of various methods and limiters on the same problem.
This notebook can be found in
$CLAW/apps/notebooks/classic/acoustics_1d_example1/acoustics_1d_example1.ipynb
The 1-dimensional acoustics equations \(q_t + Aq_x = 0\) is solved in the interval \(-1 \leq x \leq 1\) with extrapolation boundary conditions.
The density rho and bulk modulus K are specified in
setrun.py but can be changed below.
Have plots appear inline in notebook:
In [1]:
%pylab inline
Populating the interactive namespace from numpy and matplotlib
Check that the CLAW environment variable is set. (It must be set in the Unix shell before starting the notebook server).
In [2]:
import os
try:
CLAW = os.environ['CLAW']
print "Using Clawpack from ", CLAW
except:
print "*** Environment variable CLAW must be set to run code"
Using Clawpack from /Users/rjl/git/clawpack
Module with functions used to execute system commands and capture output:¶
In [3]:
from clawpack.clawutil import nbtools
Compile the code:¶
In [4]:
nbtools.make_exe(new=True) # new=True ==> force recompilation of all code
Executing shell command: make new
Done... Check this file to see output:
Make documentation files:¶
In [5]:
nbtools.make_htmls()
See the README.html file for links to input files...
Run the code and plot results using the setrun.py and setplot.py files in this directory:¶
First create data files needed for the Fortran code, using parameters specified in setrun.py:
In [6]:
nbtools.make_data(verbose=False)
Now run the code and produce plots. Specifying a label insures the resulting plot directory will persist after later runs are done below.
In [7]:
outdir,plotdir = nbtools.make_output_and_plots(label='1')
Executing shell command: make output OUTDIR=_output_1
Done... Check this file to see output:
Executing shell command: make plots OUTDIR=_output_1 PLOTDIR=_plots_1
Done... Check this file to see output:
View plots created at this link:
Display the animation inline:¶
Clicking on the _PlotIndex link above, you can view an animation of the results.
After creating all png files in the _plots directory, these can also be combined in an animation that is displayed inline:
In [8]:
import clawpack.visclaw.JSAnimation.JSAnimation_frametools as J
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim
Out[8]:
Adjust some parameters to explore the methods in Clawpack¶
The animation above was computed using the default parameter values
specified in setrun.py, which specified using the high-resolution
method with the MC limiter. See the README.html file for a link to
setrun.py.
We can adjust the parameters by reading in the default values, changing one or more, and then writing the data out for the Fortran code to use:
In [9]:
import setrun
rundata = setrun.setrun()
print "The left boundary condition is currently set to ",rundata.clawdata.bc_lower
print "The right boundary condition is currently set to ",rundata.clawdata.bc_upper
The left boundary condition is currently set to ['extrap']
The right boundary condition is currently set to ['extrap']
Periodic boundary conditions¶
Change the boundary conditions and write out the data. Then rerun the code.
In [11]:
rundata.clawdata.bc_lower = ['periodic']
rundata.clawdata.bc_upper = ['periodic']
# also extend the time over which the solution is computed and plotted:
rundata.clawdata.num_output_times = 40
rundata.clawdata.tfinal = 2.0
rundata.write()
outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim
Out[11]:
Reflecting wall boundary conditions¶
Change the boundary conditions and write out the data. Then rerun the code.
In [12]:
rundata.clawdata.bc_lower = ['wall']
rundata.clawdata.bc_upper = ['wall']
rundata.write()
outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim
Out[12]:
Things to try:¶
- Use
print rundata.clawdatato see what parameters can be set. Also printrundata.probdatato see what problem-specific paramaters are available. - Change the density
rhoor the bulk modulusKand see what effect this has.
In [13]:
def print_material():
rho = rundata.probdata.rho
K = rundata.probdata.K
Z = sqrt(K*rho)
c = sqrt(K/rho)
print "The density rho = %g and bulk modulus %g give" % (rho,K)
print " speed of sound c = %g" % c
print " impedance Z = %g" % Z
print_material()
The density rho = 1 and bulk modulus 4 give
speed of sound c = 2
impedance Z = 2
In [14]:
rundata.probdata.rho = 4.
rundata.write()
print_material()
outdir, plotdir = nbtools.make_output_and_plots(verbose=False)
anim = J.make_anim(plotdir, figno=1, figsize=(6,4))
anim
The density rho = 4 and bulk modulus 4 give
speed of sound c = 1
impedance Z = 4
Out[14]: